Libraries

   if(any(grepl("package:dvutils", search()))) detach("package:dvutils", unload=TRUE)
   library(dvutils)
   library(stringr)
   library(ggplot2)
   library(dplyr)
   library(tidyr)
   library(tibble)
   library(readr)
   library(maps)
   library(glue)
   library(readxl)
   library(rlang)
   
   library(phyloseq)

   library(knitr)
   library(rmdformats)

   library(patchwork)

   # library(pr2database)
   # data("pr2")

Markdown

Read the colors

class_colors <- read_excel("colors.xlsx", sheet = "classes")

class_colors <- structure(class_colors$color_name, .Names = class_colors$class)

genus_colors <- read_excel("colors.xlsx", sheet = "genus")

genus_colors <- structure(genus_colors$color_name, .Names = genus_colors$genus)

Initialize constants

# rm(list = ls())
sample_type <- c("18S filter", "16S plastid", "18S sort")
taxon_selected <- c("Chlorophyta", "Cryptophyta", "Rhodophyta", "Haptophyta", "Ochrophyta")

ps <- list()
long <- list()

Read the data from the database (only done once) - DO NOT RUN !

temp <- list()

# Export as specific data set as a phyloseq file

for (i in 1:3) {
    
    temp <- metapr2_export_asv(dataset_id_selected = i + 15, taxo_level = division, taxo_name = taxon_selected, export_phyloseq = TRUE, export_long_xls = FALSE, 
        export_wide_xls = TRUE, export_fasta = TRUE, export_sample_xls = TRUE, boot_min = 90, boot_level = supergroup_boot, directory = "../dada2/")
    
    
}


rm(temp)

Read the Phyloseq files from the disk

• Only keep station 6
• Do not consider TFF samples

for (i in 1:3) {
    
    ps[[sample_type[i]]] <- readRDS(str_c("../dada2/phyloseq_metapr2_asv_set_", i + 15, "_photo.rda"))
    
    
    # Filter out station 14
    ps[[sample_type[i]]] <- subset_samples(ps[[sample_type[i]]], station_id == "6")
    
    # Filter out TFF concentrated samples
    if (i == 3) {
        ps[[sample_type[i]]] <- subset_samples(ps[[sample_type[i]]], is.na(sample_concentration))
    }
    
    
    ps[[sample_type[i]]] <- ps[[sample_type[i]]] %>% # Filter ASVs that are absent
    filter_taxa(function(x) sum(x) > 0, TRUE) %>% # Filter Paraphysomonas
    subset_taxa(family != "Chrysophyceae_Clade-F")
    
    # Label samples with date
    sample_data(ps[[sample_type[i]]])$sample_label <- as.character(sample_data(ps[[sample_type[i]]])$date)
    
    # Add TFF to sample name sample_data(ps[[one_sample_type]])$sample_label <- str_c(sample_data(ps[[one_sample_type]])$metadata_code,
    # str_replace_na(sample_data(ps[[one_sample_type]])$sample_concentration, ''), sep = '.' )
    
    print(glue("Phyloseq - {sample_type[i]}"))
    print(ps[[sample_type[i]]])
    print(glue("Mean number of reads: {mean(sample_sums(ps[[ sample_type[i] ]]))}"))
    cat("============================\n")
    mean(sample_sums(ps[[sample_type[i]]]))
    
}

Phyloseq - 18S filter
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 681 taxa and 120 samples ]
sample_data() Sample Data:       [ 120 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 681 taxa by 8 taxonomic ranks ]
Mean number of reads: 27029.1083333333
============================
Phyloseq - 16S plastid
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 417 taxa and 100 samples ]
sample_data() Sample Data:       [ 100 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 417 taxa by 8 taxonomic ranks ]
Mean number of reads: 35058.32
============================
Phyloseq - 18S sort
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 195 taxa and 40 samples ]
sample_data() Sample Data:       [ 40 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 195 taxa by 8 taxonomic ranks ]
Mean number of reads: 31960.325
============================

sample_variables(ps[["18S filter"]])

 [1] "sample_id"                   "file_name"                   "NCBI_run"                    "sample_name"                
 [5] "sample_code"                 "metadata_code"               "DNA_RNA"                     "fraction_name"              
 [9] "fraction_name_original"      "fraction_min"                "fraction_max"                "reads_total"                
[13] "metadata_id"                 "project"                     "station_id"                  "station_id_num"             
[17] "year"                        "date"                        "time"                        "season"                     
[21] "depth_level"                 "depth"                       "substrate"                   "latitude"                   
[25] "longitude"                   "site_name"                   "country"                     "oceanic_region"             
[29] "temperature"                 "salinity"                    "pH"                          "O2"                         
[33] "fluorescence"                "Chla"                        "Chla_0.2_3.um"               "NO2"                        
[37] "NO3"                         "PO4"                         "Si"                          "bact_ml"                    
[41] "peuk_ml"                     "neuk_ml"                     "crypto_ml"                   "dataset_code"               
[45] "dataset_name"                "processing_pipeline_metapr2" "processing_date"             "sequencing_technology"      
[49] "sequencing_type"             "sequencing_company"          "region"                      "ecosystem"                  
[53] "lat_min"                     "lat_max"                     "long_min"                    "long_max"                   
[57] "data_available"              "bioproject_accession"        "bioproject_link"             "sample_label"               

mean(sample_sums(ps[["18S filter"]]))

[1] 27029.11

Normalize and transform to long form

for (one_sample_type in sample_type_surface) {
    
    cat(one_sample_type)
    cat("\n============================")
    
    ps[[one_sample_type]] <- phyloseq_normalize_median(ps[[one_sample_type]])
    
    long[[one_sample_type]] <- phyloseq_transform_to_long(ps[[one_sample_type]])
    
    cat("\n")
    
}

18S filter surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 383 taxa and 50 samples ]
sample_data() Sample Data:       [ 50 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 383 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  23194
16S plastid surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 236 taxa and 40 samples ]
sample_data() Sample Data:       [ 40 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 236 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  28885
18S sort surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 114 taxa and 16 samples ]
sample_data() Sample Data:       [ 16 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 114 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  27275
18S filter 0.2 um surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 187 taxa and 17 samples ]
sample_data() Sample Data:       [ 17 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 187 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  14645
18S filter 3 um surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 225 taxa and 18 samples ]
sample_data() Sample Data:       [ 18 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 225 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  24783
18S filter 20 um surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 199 taxa and 15 samples ]
sample_data() Sample Data:       [ 15 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 199 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  29427
16S plastid 0.2 um surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 112 taxa and 11 samples ]
sample_data() Sample Data:       [ 11 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 112 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  28471
16S plastid 3 um surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 142 taxa and 16 samples ]
sample_data() Sample Data:       [ 16 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 142 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  31801
16S plastid 20 um surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 120 taxa and 13 samples ]
sample_data() Sample Data:       [ 13 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 120 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  25641
18S sort pico surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 65 taxa and 8 samples ]
sample_data() Sample Data:       [ 8 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 65 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  33258
18S sort nano surface
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 70 taxa and 8 samples ]
sample_data() Sample Data:       [ 8 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 70 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  21650

List of classes

classes <- data.frame(division = NULL, class = NULL)

for (one_sample_type in sample_type_surface[1:3]) {
    classes <- bind_rows(classes, select(long[[one_sample_type]], division, class))
}

classes <- classes %>% distinct(division, class) %>% arrange(division, class)

# print(classes, n=50)
knitr::kable(classes)

Treemaps

treemap_class <- list()
treemap_genus <- list()

for (one_sample_type in sample_type_surface) {
    
    treemap_class[[one_sample_type]] <- phyloseq_long_treemap(long[[one_sample_type]], division, class, str_c(one_sample_type, " - Class"), colors = class_colors)
    treemap_genus[[one_sample_type]] <- phyloseq_long_treemap(long[[one_sample_type]], class, genus, str_c(one_sample_type, " - Genus"), colors = genus_colors)
    
}

ASVs bargraphs

bargraph_asv <- list()
bargraph_species <- list()

for (one_sample_type in sample_type_surface) {
    
    bargraph_asv[[one_sample_type]] <- phyloseq_long_bargraph(long[[one_sample_type]], text_scaling = 0.75, n_bars = 20, use_asv = TRUE, title = one_sample_type, 
        taxo_level_fill = class, taxo_colors_fill = class_colors)
    bargraph_species[[one_sample_type]] <- phyloseq_long_bargraph(long[[one_sample_type]], text_scaling = 0.75, n_bars = 20, use_asv = FALSE, title = one_sample_type, 
        taxo_level_fill = class, taxo_colors_fill = class_colors)
    
}

## Barplots per sample

bargraph_sample <- list()

for (one_sample_type in sample_type_surface[4:11]) {
    
    bargraph_sample[[one_sample_type]] <- plot_bar(ps[[one_sample_type]], x = "sample_label", fill = "class") + geom_bar(aes(color = class, fill = class), 
        stat = "identity", position = "stack") + ggtitle(str_c("Class level - ", one_sample_type)) + theme(axis.text.y = element_text(size = 10)) + 
        theme(axis.text.x = element_text(angle = 0, hjust = 0.5)) + coord_flip() + scale_fill_viridis_d() + scale_color_viridis_d()
    print(bargraph_sample[[one_sample_type]])
}

Heatmaps

Filter : Dictyo / Chrysophyceae /Pelago /Bacili / Crypto / Pyrami / Mamiello

heatmap_class_abundant <- list()

for (one_sample_type in sample_type_surface[4:11]) {
    
    ps_heat <- tax_glom(ps[[one_sample_type]], taxrank = "class") %>% phyloseq_filter_abundant_taxa(fraction_min = 0.1)
    
    heatmap_class_abundant[[one_sample_type]] <- plot_heatmap(ps_heat, method = "NMDS", distance = "bray", taxa.label = "class", taxa.order = "division", 
        sample.label = "sample_label", sample.order = "sample_label", low = "beige", high = "red", na.value = "beige", trans = NULL, title = one_sample_type) + 
        xlab("") + ylab("")
    
    # plot(heatmap(otu_table(ps_heat)))
    
    print(heatmap_class_abundant[[one_sample_type]])
    
}

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 7 taxa and 17 samples ]
sample_data() Sample Data:       [ 17 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 7 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  13806

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 7 taxa and 18 samples ]
sample_data() Sample Data:       [ 18 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 7 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  23366

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 5 taxa and 15 samples ]
sample_data() Sample Data:       [ 15 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 5 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  27860

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 4 taxa and 11 samples ]
sample_data() Sample Data:       [ 11 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 4 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  27592

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 5 taxa and 16 samples ]
sample_data() Sample Data:       [ 16 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 5 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  31460

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 4 taxa and 13 samples ]
sample_data() Sample Data:       [ 13 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 4 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  24881

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 4 taxa and 8 samples ]
sample_data() Sample Data:       [ 8 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 4 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  32016

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 4 taxa and 8 samples ]
sample_data() Sample Data:       [ 8 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 4 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  21054

class_selected <- c("Dictyochophyceae", "Chrysophyceae", "Pelagophyceae", "Bacillariophyta", "Prymnesiophyceae", "Cryptophyceae", "Pyramimonadophyceae", 
    "Mamiellophyceae")

heatmap_class_selected <- list()

for (one_sample_type in sample_type_surface[4:11]) {
    # for (one_sample_type in sample_type_surface[4:4]) {
    reads_max <- case_when(str_detect(one_sample_type, "18S filter") ~ 10000, str_detect(one_sample_type, "16S plastid") ~ 30000, str_detect(one_sample_type, 
        "18S sort") ~ 10000)
    
    if (str_detect(one_sample_type, "18S sort")) {
        class_selected <- c("Pelagophyceae", "Bacillariophyta", "Prymnesiophyceae", "Cryptophyceae", "Mamiellophyceae")
    }
    
    ps_heat <- tax_glom(ps[[one_sample_type]], taxrank = "class") %>% subset_taxa(class %in% class_selected)
    
    # Try to order by division and class...  tax_table <- data.frame(tax_table(ps_heat)@.Data) taxa_names(ps_heat) <- str_c(tax_table$division,
    # tax_table$class, sep=' - ')
    
    gg <- plot_heatmap(ps_heat, method = "NMDS", distance = "bray", taxa.label = "class", taxa.order = "class", sample.label = "sample_label", sample.order = "sample_label", 
        low = "beige", high = "red", na.value = "gray95", trans = NULL, title = one_sample_type) + theme(axis.text.x = element_text(angle = 45, 
        hjust = 1, vjust = 1), axis.title.x = element_blank(), axis.title.y = element_blank(), panel.background = element_rect(fill = "grey50"), 
        panel.grid.major = element_blank()) + scale_fill_gradient(limits = c(0, reads_max), low = "beige", high = "red", na.value = "red")
    
    if (str_detect(one_sample_type, "18S filter 0.2")) 
        gg <- gg + expand_limits(y = c(0.5, 9)) + geom_vline(xintercept = c(3.5, 9.5, 12.5)) + annotate(geom = "text", label = c("summer", "fall", 
            "spring", "summer"), x = c(2, 6.5, 11, 15), y = 8.65, color = "white", fill = "grey50", hjust = 0.5, vjust = 0, size = 4.5, label.r = unit(0.2, 
            "lines"))
    
    
    if (one_sample_type == "18S filter 3 um surface") 
        gg <- gg + expand_limits(y = c(0.5, 9)) + geom_vline(xintercept = c(4.5, 10.5, 13.5)) + annotate(geom = "text", label = c("summer", "fall", 
            "spring", "summer"), x = c(2.5, 7.5, 12, 16), y = 8.65, color = "white", fill = "grey50", hjust = 0.5, vjust = 0, size = 4.5, label.r = unit(0.2, 
            "lines"))
    
    if (one_sample_type == "18S filter 20 um surface") 
        gg <- gg + expand_limits(y = c(0.5, 9)) + geom_vline(xintercept = c(4.5, 9.5)) + annotate(geom = "text", label = c("summer", "fall", "summer"), 
            x = c(2.5, 7, 12.5), y = 8.65, color = "white", fill = "grey50", hjust = 0.5, vjust = 0, size = 4.5, label.r = unit(0.2, "lines"))
    
    
    # plot(heatmap(otu_table(ps_heat))) , trans = scales::log_trans(10)
    
    print(gg)
    
    heatmap_class_selected[[one_sample_type]] <- gg
    
}

heatmap_genus_abundant <- list()

for (one_sample_type in sample_type_surface[4:11]) {
    
    ps_heat <- tax_glom(ps[[one_sample_type]], taxrank = "species") %>% phyloseq_filter_abundant_taxa(fraction_min = 0.1)
    
    heatmap_genus_abundant[[one_sample_type]] <- plot_heatmap(ps_heat, method = "NMDS", distance = "bray", taxa.label = "species", taxa.order = "division", 
        sample.label = "sample_label", sample.order = "sample_label", low = "beige", high = "red", na.value = "beige", title = one_sample_type)
    
    print(heatmap_genus_abundant[[one_sample_type]])
    
}

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 10 taxa and 17 samples ]
sample_data() Sample Data:       [ 17 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 10 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  11702

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 8 taxa and 18 samples ]
sample_data() Sample Data:       [ 18 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 8 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  20199

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 11 taxa and 15 samples ]
sample_data() Sample Data:       [ 15 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 11 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  26302

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 6 taxa and 11 samples ]
sample_data() Sample Data:       [ 11 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 6 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  27286

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 9 taxa and 16 samples ]
sample_data() Sample Data:       [ 16 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 9 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  30697

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 8 taxa and 13 samples ]
sample_data() Sample Data:       [ 13 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 8 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  24093

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 5 taxa and 8 samples ]
sample_data() Sample Data:       [ 8 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 5 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  28373

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 6 taxa and 8 samples ]
sample_data() Sample Data:       [ 8 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 6 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  19513

species_selected <- c("Thalassiosira_minima", "Fragilariopsis_cylindrus", "Minidiscus_sp.", "Chaetoceros_neogracilis", "Porosira_glacialis", "Coretrhon_inerme", 
    "Palmaria_palmata", "Pseudo-nitzchia_seriata", "Micromonas_polaris", "Micromonas_clade_B3", "Bathycoccus_prasinos", "Pyramimonas_gelidicola", 
    "Geminigera_cryophila", "Pelagophyceae_XXX_sp.", "Phaeocystis_antarctica")


heatmap_species_selected <- list()

for (one_sample_type in sample_type_surface[4:11]) {
    # for (one_sample_type in sample_type_surface[4:4]) {
    reads_max <- case_when(str_detect(one_sample_type, "18S filter") ~ 10000, str_detect(one_sample_type, "16S plastid") ~ 5000, str_detect(one_sample_type, 
        "18S sort") ~ 10000)
    if (str_detect(one_sample_type, "18S sort")) {
        species_selected <- c("Thalassiosira_minima", "Fragilariopsis_cylindrus", "Minidiscus_sp.", "Chaetoceros_neogracilis", "Porosira_glacialis", 
            "Coretrhon_inerme", "Palmaria_palmata", "Pseudo-nitzchia_seriata", "Micromonas_polaris", "Geminigera_cryophila", "Pelagophyceae_XXX_sp.", 
            "Phaeocystis_antarctica")
    }
    
    
    
    ps_heat <- tax_glom(ps[[one_sample_type]], taxrank = "species") %>% subset_taxa(species %in% species_selected)
    
    # Try to order by division and species...  tax_table <- data.frame(tax_table(ps_heat)@.Data) taxa_names(ps_heat) <- str_c(tax_table$division,
    # tax_table$species, sep=' - ')
    
    gg <- plot_heatmap(ps_heat, method = "NMDS", distance = "bray", taxa.label = "species", taxa.order = "class", sample.label = "sample_label", 
        sample.order = "sample_label", low = "beige", high = "red", na.value = "gray95", trans = NULL, title = one_sample_type) + theme(axis.text.x = element_text(angle = 45, 
        hjust = 1, vjust = 1), axis.title.x = element_blank(), axis.title.y = element_blank(), panel.background = element_rect(fill = "grey50"), 
        panel.grid.major = element_blank()) + scale_fill_gradient(limits = c(0, reads_max), low = "beige", high = "red", na.value = "red")
    
    if (str_detect(one_sample_type, "18S filter 0.2")) 
        gg <- gg + geom_vline(xintercept = c(3.5, 9.5, 12.5))
    
    if (one_sample_type == "18S filter 3 um surface") 
        gg <- gg + geom_vline(xintercept = c(4.5, 10.5, 13.5))
    
    if (one_sample_type == "18S filter 20 um surface") 
        gg <- gg + geom_vline(xintercept = c(4.5, 9.5))
    
    
    
    # plot(heatmap(otu_table(ps_heat))) , trans = scales::log_trans(10)
    
    print(gg)
    
    heatmap_species_selected[[one_sample_type]] <- gg
    
}

Apply to all samples

nmds_sample <- list()

for (one_sample_type in sample_type_surface[1:3]) {
    
    one_ps <- ps[[one_sample_type]] %>% tax_glom(taxrank = "species")
    
    nmds_sample[[one_sample_type]] <- phyloseq_nmds(one_ps, title = one_sample_type, sample_color = season, sample_shape = fraction_name_original, 
        taxo_level = class, taxo_colors = class_colors)
}

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1714216 
Run 1 stress 0.1758561 
Run 2 stress 0.1717706 
... Procrustes: rmse 0.007217713  max resid 0.02514368 
Run 3 stress 0.1759724 
Run 4 stress 0.1714241 
... Procrustes: rmse 0.0003112062  max resid 0.001279752 
... Similar to previous best
Run 5 stress 0.1713535 
... New best solution
... Procrustes: rmse 0.003791301  max resid 0.01300957 
Run 6 stress 0.1713543 
... Procrustes: rmse 0.0002738691  max resid 0.00105073 
... Similar to previous best
Run 7 stress 0.1714036 
... Procrustes: rmse 0.003856787  max resid 0.01498517 
Run 8 stress 0.1758565 
Run 9 stress 0.1760002 
Run 10 stress 0.1713761 
... Procrustes: rmse 0.002503714  max resid 0.01167642 
Run 11 stress 0.190753 
Run 12 stress 0.1758566 
Run 13 stress 0.1778862 
Run 14 stress 0.1776749 
Run 15 stress 0.1937439 
Run 16 stress 0.1713777 
... Procrustes: rmse 0.002296881  max resid 0.01177999 
Run 17 stress 0.1759718 
Run 18 stress 0.1713972 
... Procrustes: rmse 0.002984647  max resid 0.0119656 
Run 19 stress 0.1866215 
Run 20 stress 0.1904543 
*** Solution reached
[1] 3

Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance) 

global Multidimensional Scaling using monoMDS

Data:     wisconsin(sqrt(veganifyOTU(physeq))) 
Distance: bray 

Dimensions: 2 
Stress:     0.1713535 
Stress type 1, weak ties
Two convergent solutions found after 20 tries
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))' 

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1477193 
Run 1 stress 0.151694 
Run 2 stress 0.1515777 
Run 3 stress 0.1720669 
Run 4 stress 0.1476132 
... New best solution
... Procrustes: rmse 0.008944154  max resid 0.03744427 
Run 5 stress 0.1743639 
Run 6 stress 0.1651081 
Run 7 stress 0.1531319 
Run 8 stress 0.1490493 
Run 9 stress 0.165109 
Run 10 stress 0.1651081 
Run 11 stress 0.1925171 
Run 12 stress 0.162316 
Run 13 stress 0.1737626 
Run 14 stress 0.1720647 
Run 15 stress 0.1519371 
Run 16 stress 0.1736789 
Run 17 stress 0.1720671 
Run 18 stress 0.1874055 
Run 19 stress 0.150668 
Run 20 stress 0.1742395 
*** No convergence -- monoMDS stopping criteria:
    20: stress ratio > sratmax
[1] 3

Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance) 

global Multidimensional Scaling using monoMDS

Data:     wisconsin(sqrt(veganifyOTU(physeq))) 
Distance: bray 

Dimensions: 2 
Stress:     0.1476132 
Stress type 1, weak ties
No convergent solutions - best solution after 20 tries
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))' 

Square root transformation
Wisconsin double standardization
Run 0 stress 0.05090316 
Run 1 stress 0.0524764 
Run 2 stress 0.05145132 
Run 3 stress 0.05736205 
Run 4 stress 0.05466462 
Run 5 stress 0.05537261 
Run 6 stress 0.06026908 
Run 7 stress 0.0524764 
Run 8 stress 0.05161248 
Run 9 stress 0.05185409 
Run 10 stress 0.05723175 
Run 11 stress 0.05145131 
Run 12 stress 0.05945005 
Run 13 stress 0.05440615 
Run 14 stress 0.05202902 
Run 15 stress 0.3065991 
Run 16 stress 0.06090273 
Run 17 stress 0.05950989 
Run 18 stress 0.05188958 
Run 19 stress 0.05220003 
Run 20 stress 0.05185435 
*** No convergence -- monoMDS stopping criteria:
    20: stress ratio > sratmax
[1] 3

Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance) 

global Multidimensional Scaling using monoMDS

Data:     wisconsin(sqrt(veganifyOTU(physeq))) 
Distance: bray 

Dimensions: 2 
Stress:     0.05090316 
Stress type 1, weak ties
No convergent solutions - best solution after 20 tries
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))' 

Filter the data

sample_type_profile <- str_c(sample_type_all, " profile")

for (one_sample_type in sample_type_all) {
    
    one_sample_type_profile = str_c(one_sample_type, " profile")
    
    ps[[one_sample_type_profile]] <- ps[[one_sample_type]] %>% subset_samples(station_id == "6") %>% subset_samples(date == "2015-01-16") %>% filter_taxa(function(x) sum(x) > 
        0, TRUE) %>% phyloseq_normalize_median()
    
    # Rename the samples
    
    if (!str_detect(one_sample_type_profile, "sort")) {
        sample_names(ps[[one_sample_type_profile]]) <- str_c(as.character(sample_data(ps[[one_sample_type_profile]])$depth), "m", as.character(sample_data(ps[[one_sample_type_profile]])$fraction_min), 
            sep = "_")
    } else {
        sample_names(ps[[one_sample_type_profile]]) <- str_c(as.character(sample_data(ps[[one_sample_type_profile]])$depth), "m", as.character(sample_data(ps[[one_sample_type_profile]])$fraction_name), 
            sep = "_")
    }
    
    
    print(glue("Phyloseq - {one_sample_type_profile}"))
    print(ps[[one_sample_type_profile]])
    cat("============================\n")
    
}

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 172 taxa and 15 samples ]
sample_data() Sample Data:       [ 15 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 172 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  43461Phyloseq - 18S filter profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 172 taxa and 15 samples ]
sample_data() Sample Data:       [ 15 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 172 taxa by 8 taxonomic ranks ]
============================

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 107 taxa and 15 samples ]
sample_data() Sample Data:       [ 15 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 107 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  41130Phyloseq - 16S plastid profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 107 taxa and 15 samples ]
sample_data() Sample Data:       [ 15 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 107 taxa by 8 taxonomic ranks ]
============================

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 94 taxa and 10 samples ]
sample_data() Sample Data:       [ 10 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 94 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  29304Phyloseq - 18S sort profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 94 taxa and 10 samples ]
sample_data() Sample Data:       [ 10 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 94 taxa by 8 taxonomic ranks ]
============================

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 90 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 90 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  30257Phyloseq - 18S filter 0.2 um profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 90 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 90 taxa by 8 taxonomic ranks ]
============================

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 84 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 84 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  31205Phyloseq - 18S filter 3 um profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 84 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 84 taxa by 8 taxonomic ranks ]
============================

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 121 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 121 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  58797Phyloseq - 18S filter 20 um profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 121 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 121 taxa by 8 taxonomic ranks ]
============================

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 46 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 46 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  40267Phyloseq - 16S plastid 0.2 um profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 46 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 46 taxa by 8 taxonomic ranks ]
============================

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 62 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 62 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  46172Phyloseq - 16S plastid 3 um profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 62 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 62 taxa by 8 taxonomic ranks ]
============================

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 74 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 74 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  47020Phyloseq - 16S plastid 20 um profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 74 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 60 sample variables ]
tax_table()   Taxonomy Table:    [ 74 taxa by 8 taxonomic ranks ]
============================

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 63 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 63 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  29412Phyloseq - 18S sort pico profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 63 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 63 taxa by 8 taxonomic ranks ]
============================

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 46 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 46 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  28502Phyloseq - 18S sort nano profile
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 46 taxa and 5 samples ]
sample_data() Sample Data:       [ 5 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 46 taxa by 8 taxonomic ranks ]
============================

Barplots per depth

bargraph_sample <- list()

for (one_sample_type in sample_type_profile[4:11]) {
    
    bargraph_sample[[one_sample_type]] <- plot_bar(ps[[one_sample_type]], x = "depth", fill = "class") + geom_bar(aes(color = class, fill = class), 
        stat = "identity", position = "stack") + ggtitle(str_c("Class level - ", one_sample_type)) + theme(axis.text.y = element_text(size = 10)) + 
        theme(axis.text.x = element_text(angle = 0, hjust = 0.5)) + coord_flip() + # scale_fill_viridis_d() + scale_color_viridis_d() +
    scale_fill_manual(values = class_colors) + scale_color_manual(values = class_colors) + scale_x_reverse() + xlab("Depth (m)") + ylab("Reads") + 
        theme_bw()
    
    print(bargraph_sample[[one_sample_type]])
}

Sample by sample comparison

• Each dot corresponds to the relative contribution of the taxonomic level considered for a sample for which both 18S and 16S has been performed

long_species_sample <- list()
long_class_sample <- list()

for (one_sample_type in sample_type) {
    ps_one <- phyloseq_normalize_percent(ps[[one_sample_type]])
    
    long_one <- phyloseq_transform_to_long(ps_one) %>% filter(!is.na(fraction_name_original)) %>% group_by(date, station_id, depth, fraction_name_original, 
        division, class, order, family, genus, species) %>% summarize(n_reads = sum(n_reads)) %>% ungroup()
    
    long_species_sample[[one_sample_type]] <- long_one
    
    long_one <- phyloseq_transform_to_long(ps_one) %>% filter(!is.na(fraction_name_original)) %>% group_by(date, station_id, depth, fraction_name_original, 
        division, class) %>% summarize(n_reads = sum(n_reads)) %>% ungroup()
    
    long_class_sample[[one_sample_type]] <- long_one
}

long_compare_species_18S_16S <- full_join(rename(long_species_sample[["18S filter"]], n_reads_18S_filter = n_reads), rename(long_species_sample[["16S plastid"]], 
    n_reads_16S_plastid = n_reads))

ggplot(long_compare_species_18S_16S, aes(x = n_reads_18S_filter, y = n_reads_16S_plastid)) + geom_point(aes(color = class)) + scale_color_manual(values = class_colors) + 
    ggtitle("Species")

long_compare_class_18S_16S <- long_compare_species_18S_16S %>% group_by(date, station_id, depth, fraction_name_original, division, class) %>% summarize(n_reads_18S_filter = sum(n_reads_18S_filter, 
    na.rm = TRUE), n_reads_16S_plastid = sum(n_reads_16S_plastid, na.rm = TRUE))

ggplot(long_compare_class_18S_16S, aes(x = n_reads_18S_filter, y = n_reads_16S_plastid)) + geom_point(aes(color = class)) + scale_color_manual(values = class_colors) + 
    ggtitle("Classes")

Define upset function

• Needs only to have 0 and 1 in the matrix

plot_upset_long <- function(long_comparison, file_name) {
    
    long_upset <- long_comparison %>% mutate_if(is.numeric, funs(case_when(is.na(.) ~ 0, TRUE ~ 1))) %>% as.data.frame()
    
    fig_upset <- UpSetR::upset(long_upset, empty.intersections = "on", number.angles = 30, point.size = 12, line.size = 4, mainbar.y.label = "Taxa number", 
        sets.x.label = "Taxa #", order.by = "freq", text.scale = c(5), mb.ratio = c(0.7, 0.3))
    
    # text.scale = c(intersection size title, intersection size tick labels, set size title, set size tick labels, set names, numbers above bars)
    
    print(fig_upset)
    png(file_name, width = 2000, height = 1200, units = "px")
    print(fig_upset)
    dev.off()
    
    return(fig_upset)
}

Global comparison of species found with the 3 methods

• Compute table of number of samples for each species (rows) vs the three methods (columns)

• Only keep

  • 2015 samples because it is the only dataset for which we have the three types of samples
  • 0.2 and 3 um fractions (to be comparable with sorting)

• Remove species that contains _X

long_species <- list()
long_genus <- list()
long_class <- list()
fig_upset_list <- list()

for (one_sample_type in sample_type) {
    long_one <- long_species_sample[[one_sample_type]] %>% filter(date >= as.Date("2015-01-01")) %>% filter(!str_detect(species, "_X")) %>% filter(fraction_name_original != 
        "> 20 um") %>% group_by(division, class, order, family, genus, species) %>% summarize(n_reads = sum(n_reads), n_samples = n()) %>% rename_at(vars(contains("n_")), 
        funs(str_c(., str_replace_all(one_sample_type, " ", "_"), sep = "_"))) %>% ungroup()
    
    long_species[[one_sample_type]] <- long_one
    
    long_one <- long_species_sample[[one_sample_type]] %>% filter(date >= as.Date("2015-01-01")) %>% filter(!str_detect(species, "_X")) %>% filter(fraction_name_original != 
        "> 20 um") %>% group_by(division, class, order, family, genus) %>% summarize(n_reads = sum(n_reads), n_samples = n()) %>% rename_at(vars(contains("n_")), 
        funs(str_c(., str_replace_all(one_sample_type, " ", "_"), sep = "_"))) %>% ungroup()
    
    long_genus[[one_sample_type]] <- long_one
    
    long_one <- long_class_sample[[one_sample_type]] %>% filter(date >= as.Date("2015-01-01")) %>% filter(fraction_name_original != "> 20 um") %>% 
        group_by(division, class) %>% summarize(n_reads = sum(n_reads), n_samples = n()) %>% rename_at(vars(contains("n_")), funs(str_c(., str_replace_all(one_sample_type, 
        " ", "_"), sep = "_"))) %>% ungroup()
    
    long_class[[one_sample_type]] <- long_one
    
}

long_class_comparison <- long_class %>% purrr::reduce(full_join)

long_genus_comparison <- long_genus %>% purrr::reduce(full_join)

long_species_comparison <- long_species %>% purrr::reduce(full_join)

# Class ----------

cat("Class not found in one type of sample\n")

Class not found in one type of sample

cat("18S filter\n")

18S filter

kable(filter(long_class_comparison, is.na(n_samples_18S_filter)))

cat("18S sort\n")

18S sort

kable(filter(long_class_comparison, is.na(n_samples_18S_sort)))

cat("16S filter\n")

16S filter

kable(filter(long_class_comparison, is.na(n_samples_16S_plastid)))

cat("Genus  found in the three types of samples\n")

Genus  found in the three types of samples

cat("18S filter\n")

18S filter

kable(filter(long_genus_comparison, (!is.na(n_samples_18S_filter) & !is.na(n_samples_18S_sort) & !is.na(n_samples_16S_plastid))))

# Genus ----------

cat("Genus only found in one type of sample\n")

Genus only found in one type of sample

cat("18S filter\n")

18S filter

kable(filter(long_genus_comparison, (is.na(n_samples_18S_sort) & is.na(n_samples_16S_plastid))))

cat("18S sort\n")

18S sort

kable(filter(long_genus_comparison, (is.na(n_samples_18S_filter) & is.na(n_samples_16S_plastid))))

cat("16S filter\n")

16S filter

kable(filter(long_genus_comparison, (is.na(n_samples_18S_filter) & is.na(n_samples_18S_sort))))

cat("Species  found in the three types of samples\n")

Species  found in the three types of samples

cat("18S filter\n")

18S filter

kable(filter(long_species_comparison, (!is.na(n_samples_18S_filter) & !is.na(n_samples_18S_sort) & !is.na(n_samples_16S_plastid))))

long_genus_comparison_upset <- long_genus_comparison %>% select(genus, contains("samples")) %>% rename_at(vars(contains("n_")), funs(str_replace_all(., 
    c(n_samples_ = "", `_` = " ", plastid = "filter"))))

fig_upset_list[["method genus"]] <- plot_upset_long(long_genus_comparison_upset, "../fig/upset_method_genus.png")

# Species ----------


cat("Species only found in one type of sample\n")

Species only found in one type of sample

cat("18S filter\n")

18S filter

kable(filter(long_species_comparison, (is.na(n_samples_18S_sort) & is.na(n_samples_16S_plastid))))

cat("18S sort\n")

18S sort

kable(filter(long_species_comparison, (is.na(n_samples_18S_filter) & is.na(n_samples_16S_plastid))))

cat("16S filter\n")

16S filter

kable(filter(long_species_comparison, (is.na(n_samples_18S_filter) & is.na(n_samples_18S_sort))))

long_species_comparison_upset <- long_species_comparison %>% select(species, contains("samples")) %>% rename_at(vars(contains("n_")), funs(str_replace(., 
    "n_samples_", "")))

fig_upset_list[["method species"]] <- plot_upset_long(long_species_comparison_upset, "../fig/upset_method_species.png")

Write to Excel file

openxlsx::write.xlsx(list(long_class_comparison, long_genus_comparison, long_species_comparison), file = "../dada2/method_comparison.xlsx", sheetName = c("Class", 
    "Genus", "Species"))

Comparison of species found in the different fractions for 18S filter and 18S sorted

• Compute table of number of samples for each species (rows) vs the three fractions (columns)

• Only keep

  • 2015 samples because it is the only dataset for which we have the three types of samples
  • Genera that do not contain _X

long_species_fraction <- list()
long_class_fraction <- list()

for (one_sample_type in sample_type[c(1, 3)]) {
    long_one <- long_species_sample[[one_sample_type]] %>% filter(date >= as.Date("2015-01-01")) %>% filter(!str_detect(species, "_X")) %>% group_by(fraction_name_original, 
        division, class, order, family, genus, species) %>% summarize(n_samples = n()) %>% pivot_wider(names_from = fraction_name_original, values_from = n_samples) %>% 
        ungroup()
    
    long_species_fraction[[one_sample_type]] <- long_one
    
    long_one <- long_class_sample[[one_sample_type]] %>% filter(date >= as.Date("2015-01-01")) %>% group_by(fraction_name_original, division, class) %>% 
        summarize(n_samples = n()) %>% pivot_wider(names_from = fraction_name_original, values_from = n_samples) %>% ungroup()
    
    long_class_fraction[[one_sample_type]] <- long_one
    
}

long_class_fraction_comparison <- long_class_fraction %>% purrr::reduce(full_join)
long_species_fraction_comparison <- long_species_fraction %>% purrr::reduce(full_join)

kable(long_class_fraction_comparison)

kable(long_species_fraction_comparison)

fig_upset_list[["filter_sort species"]] <- plot_upset_long(long_species_fraction_comparison, "../fig/upset_filter_sort.png")

Comparison of species found in the different fractions for 18S filter

• Compute table of number of samples for each species (rows) vs the three fractions (columns)

• Only keep

  • Genera that do not contain _X

• We consider now ALL samples (will need to remove samples for which the > 20 um is missing)

long_species_fraction <- list()
long_genus_fraction <- list()
long_class_fraction <- list()

for (one_sample_type in sample_type[c(1)]) {
    long_one <- long_species_sample[[one_sample_type]] %>% filter(!str_detect(species, "_X")) %>% group_by(fraction_name_original, division, class, 
        order, family, genus, species) %>% summarize(n_samples = n()) %>% pivot_wider(names_from = fraction_name_original, values_from = n_samples) %>% 
        ungroup()
    
    long_species_fraction[[one_sample_type]] <- long_one
    
    long_one <- long_species_sample[[one_sample_type]] %>% filter(!str_detect(species, "_X")) %>% group_by(fraction_name_original, division, class, 
        order, family, genus) %>% summarize(n_samples = n()) %>% pivot_wider(names_from = fraction_name_original, values_from = n_samples) %>% ungroup()
    
    long_genus_fraction[[one_sample_type]] <- long_one
    
    
    long_one <- long_class_sample[[one_sample_type]] %>% group_by(fraction_name_original, division, class) %>% summarize(n_samples = n()) %>% pivot_wider(names_from = fraction_name_original, 
        values_from = n_samples) %>% ungroup()
    
    long_class_fraction[[one_sample_type]] <- long_one
    
}

long_class_fraction_comparison <- long_class_fraction[[one_sample_type]]

long_genus_fraction_comparison <- long_genus_fraction[[one_sample_type]]

long_species_fraction_comparison <- long_species_fraction[[one_sample_type]]

kable(long_class_fraction_comparison)

kable(long_genus_fraction_comparison)